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The temperature fluctuations generated by viscous dissipation in an isotropic turbulent flow are 
studied using direct numerical simulation. It is shown that their scaling with Reynolds number is at 
odds with predictions from recent investigations. The origin of the discrepancy is traced back to the 
anomalous scaling of the dissipation rate fluctuations. Phenomenological arguments are presented 
which explain the observed results. The study shows that previously proposed models underpredict 
the variance of frictional temperature fluctuations by a factor proportional to the square of the 
Taylor-scale Reynolds number. 


I. INTRODUCTION 

Two recent investigations have considered the viscous 
generation of temperature fluctuations in isotropic turbu¬ 
lent flows Hi- Even though the models fundamentally 
differ in the form of the heat production term, they pre¬ 
dict the same scaling of the wavenumber spectrum of the 
heat fluctuations in the case of a statistically steady state. 
In the second of these works concerns were expressed 
about the capability of the model to correctly predict 
the viscous heat production since it was derived using 
the Direct Interaction Approximation (DIA) Q, which 
is incapable of predicting the cumulant contributions to 
the wavenumber spectrum of the dissipation rate fluctua¬ 
tions 0 . Indeed, the temperature fluctuations generated 
by frictional heating are intimately related to the fluctua¬ 
tions of the dissipation. In the present work we will first 
revisit the derivation presented in reference 0 to pro¬ 
pose an alternative description of the physics of turbulent 
frictional heating. After that, we will check the different 
assumptions using Direct Numerical Simulations. 


II. ANALYSIS OF THE PROBLEM AND 
SCALING ARGUMENTS 


A. Governing equations 


Using the Reynolds decomposition we write the equa¬ 
tion for the fluctuation of the temperature 0 = 0 — 0, 


ae ae a^e i , 


( 1 ) 


where a bar indicates an ensemble average, Cp is the spe¬ 
cific heat, a the thermal diffusivity, and Ui an isotropic 
turbulent velocity field. The fluctuation of the dissipa¬ 
tion rate is given by 
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( 2 ) 


For notational convenience e will be used to denote the 
average dissipation rate. Since we will not use the total 
dissipation rate, the absence of a bar on this quantity 
will not introduce any ambiguity. 

The temperature distribution over different length- 
scales is given by the temperature spectrum, defined such 
that 


jEg{k)dk = 9^, (3) 

and the mean dissipation of heat fluctuations is related 
to this quantity by the relation 

J 2ak'^Eg{k)dk = eg. (4) 

The evolution-equation of Eg{k) is 

^^^=Tg{k)-Dg{k)+Pg{k), (5) 

where Tg{k) is the interscale transfer term, Dg{k) the 

temperature fluctuation dissipation term and Pg{k) the 
heat-production term. This last term is given by 

4-7rI*^ 

Pg{k) = ^{e'{k)9{-k)). (6) 

Cp 

Our interest will be, in particular, in the quantities 0^, 
eg and Eg{k). 


B. The viscous heat production as a function of 
the dissipation rate fluctuations 

In reference [^, expression m was modeled using an 
approach akin to the Direct Interaction Approximation, 
expressing {e'{k)9{—k)) as a function of velocity triad- 
interactions only. The concerns expressed in that work 
motivate us to not express this correlation as a function 
of the velocity, but to leave explicitly the dissipation rate 
fluctuations in the expression. 
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An amount of scalar 9 is advected, diffused and pro¬ 
duced. The production part, ignoring the deformation of 
the scalar blob, can be formally written according to m, 
as 

0{x,t) = 9{x,t\0) + — [ [ geix,t\y,s)e'{y, s)dyds,(7) 

where the time-integral is taken on a Lagrangian trajec¬ 
tory: 9{x,t\0) is the value of the scalar fluctuation at 
time < = 0 in the fluid particle that passes through x at 
time t. The value of 9 is thus determined by the cumula¬ 
tive heating, proportional to the dissipation rate fluctua¬ 
tion, along a trajectory. In this expression g 0 {x,t\y, s) 
is thus the Lagrangian scalar Green’s-function. The 
temperature-dissipation correlation, related to the pro¬ 
duction term in ([61), is then in a statistically steady state, 


Expression (jO)) shows that, according to our arguments, 
the viscous heat-production is proportional to the spec¬ 
trum of the dissipation rate fluctuations. It is thus this 
latter quantity that plays a major role in the determina¬ 
tion of the viscous heating. 

The spectrum E^{k) has received considerable interest 
in the past, in particular since its behavior is largely in 
disagreement with dimensional scaling within the Kol¬ 
mogorov framework Q. In particular, the latter scaling 
theory would predict the dissipation rate spectrum to 
behave as 

E^P{k) ( 12 ) 

This scaling can be obtained by assuming the velocity 
to satisfy a multi-variate Gaussian distribution and as¬ 
suming the energy spectrum to have an inertial range 
proportional to 

(13) 


{e'{x, t)9{x + r,t)) = — [ [ {gg{x + r, t\y, s)e'(a;, t)e'{y, s)) dyds.{ 8 ) 
Cp Jo J 


E{k) ~ 


We do not attempt here to derive a Lagrangian closure 
theory for this quantity, but from the above expression 
we understand that the production is related to the prod¬ 
uct of dissipation-rate fluctuations, integrated over a La¬ 
grangian correlation time. From equation ([8]), combined 
with equation (jG]), we suggest therefore that the produc¬ 
tion term can be modeled as 

Pg{k) - ^E,{k). (9) 

Cp 


Exp eriments (see @ for a review) and simulations 

show howevever that the scaling (HU) is not observed 
at high Reynolds numbers and that, instead of an increas¬ 
ing function of the wavenumber, E^(k) is a decreasing 
function of k. Yaglom proposed a model which fits the 
data qualitatively, taking into account the non-Gaussian 
character of the dissipation rate @,[11. His model pre¬ 
dicts the spectrum 

E;(2)(/c) - e^L{kL)-^+>^, (14) 


with Ee{k) defined such that 

J E,{k)dk=7^. ( 10 ) 

The timescale r(fc) in expression ([S]) represents the cor¬ 
relation time of the scalar fluctuation over a Lagrangian 
trajectory, combined with the correlation time of the 
dissipation rate fluctuations. It is at present not en¬ 
tirely clear what the functional form of this timescale is. 
The timescale of inertial range eddies of scale I ~ 1/fc, 
consistent with Kolmogorov-Obukhov scaling, is t (/ c ) ^ 
which is the Lagrangian timescale that is 
most commonly used in turbulence theory @]. However, 
for several phenomena such as Burger’s turbulence or the 
advection of a scalar in the presence of a mean scalar 
gradient, the integral timescale T (L is the 

integal lengthscale) can play a dominant role IQ. In 
the present case, as we will see, the dissipation-rate fluc¬ 
tuations are correlated at large scales, so that it is not 
excluded that the dynamics are governed by E. In order 
to leave open the possibility of both timescales to play a 
role, we introduce the hybrid timescale 

r(fc) ~ e-i/3fc-2/3(fcL)“. (II) 

The parameter a allows to consider the different possi¬ 
bilities. For a = 0 we have the Kolmogorov-Obukhov 
timescale and for a = 2/3 we find the integral timescale. 


where L is a large-scale correlation length, and fi is an 
intermittency parameter. In reference values are re¬ 
ported in the interval 0.3 < y < 0.5. Since then, a num¬ 
ber of different models have been proposed to more pre¬ 
cisely predict the physics of the dissipation rate, but these 
refinements are beyond the scope of the present investi¬ 
gation, given the accuracy obtained in the determination 
of scaling exponents at low and moderate Reynolds num¬ 
bers. 

Substituting expression (IT^ in the expression for the 
heat production term and using r(fc) ~ e ^^^k we 
find for the Gaussian estimate of the production term, 


P^^\k) 



(15) 


and this is what was obtained in reference @]. 

If we use however the scaling proposed by Yaglom m, 
combined with the timescale (HU, we find the production 
term proportional to 


^5/3u-5/3+fi+aTp+a 

Pe\k) --, (16) 

Cp 


and we recall that that y is the intermittency parameter 
in the range 0.3 < y < 0.5, and a characterizes the 
timescale (a = 0 for a Kolmogorov-Obukhov timescale, 
a = 218 for the integral timescale). 
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While getting heated on its trajectory, the fluid blob 
will also be deformed, thereby exchanging its tempera¬ 
ture variance with other lengthscales. This effect is rep¬ 
resented by the term Tg{k) in expression ([S]). We will 
investigate the scaling of the temperature spectrum, as 
in [^, by assuming the scalar transfer to be given by a 
Kovaznay-type scalar transfer model [l^, which allows 
simple analytical treatment. This model reads 

Teik) - ^ {Ee{k)E{kf/^k^/^) . (17) 

In the statistically stationary case, in the range where the 
production is important, a balance is expected between 
the transfer and the production, 


Tg{k) ^ Pe{k). (18) 


We will consider the case of unity Prandtl number, v = a, 
and we suppose the spectra to extend from ^ 1 /L to 
k^, with kri ~ ei/4j/-3/4 ^ (display 

power law behaviour throughout. The energy spectrum 
is given by m and the production spectrum is given 
either by m or m- Integrating (fT51) from k^ to k, we 
have from (II3 and 0 


Eg{k) 


£- 1 / 3 ^- 5/3 f 

JkL 


(19) 


For the Gaussian estimate of the production spec¬ 
trum m and the Kolmogorov-Obukhov time T{k) ^ 
this gives for k kh, 


E^^\k) 


£2/3^2^1/3 


( 20 ) 


For the production term taking into account intermit- 
tency effects we find 


Ef\k) 


£4/3^2/3^-573 


{{kL)^^+°‘-^/^ - l) . (21) 


larger than —1 for a = 2/3. In the former case the spec¬ 
tral flux, determined by the integral of the production 
up to wavenumber k, becomes then independent of the 
forcing exponent. For the above integrals to converge 
we need to satisfy the constraint 0 < fi < 2/3, which 
includes the experimental interval 0.3 < fj, < 0.5. 

The variance of the temperature fluctuations can be 
computed by integrating the spectra E^\k) and Eg^\k) 
between kL and kg. For the two different predictions (|^ 
and (|^ we find. 


02(4) 


eu 
72 ’ 


^(2) ^ (££)4/3 


(24) 

(25) 


and both values a = 0, 2/3 give the same scaling for the 

-( 2 ) 

variance of the temperature fluctuations 9“^ , indepen¬ 

dent of the value of /i. This is not so for the destruction 
of scalar variance. 


ee 


Pg{k)dk. 


(26) 


For the Gaussian case it gives 


^9 


. 3 / 2 , 3/2 


For the non-Gaussian estimate it yields 

^(2) ^ £"'"n q, _ q 

for a = 2/3, 


„ 5 / 3 ^ 2/3 

-72- 


65 / 3 ^ 2/3 ^ 3^/4 

cl 


(27) 


(28) 

(29) 


where Rl = UL/v and where we used e ~ U^/L with U 
the RMS velocity. The normalized variance and scalar- 
dissipation can then be defined by 


6»2 


( 1 ) 


6»2 


( 2 ) 


02 


ev 

(eL)4/3 ’ 


7(4) ^ . P 

e g3/2j^l/2 ’ 

£>2 

7(2) _ , S 


(30) 

(31) 


We now consider the different values of the exponent for 
the timescale r. We find for kL 1 


Ef\k) 


e4/3L2/3^-5/3 


for a = 0, 


( 22 ) 


independent of fx. Alternatively, it gives 


El^\k) 


e4/3L2/3fc-5/3 


(fcL))" for a = 2/3, (23) 


and the spectrum becomes thus proportional to k 
The fact that, depending on the value of a, expression 
m yields two qualitatively different spectra (one depen¬ 
dent on /t, the other not), stems from the fact that the 
forcing term (HD) is a multiscale (or fractal) forcing, with 
a power law exponent smaller than —1 for a = 0, but 


The quantities 0*2^ \ 0*2^ ^ and should become con¬ 
stant at large values of the Reynolds number, if the 
underlying assumptions are correct. The quantity ig ' 
should become constant if a = 0. However, if a = 2/3, it 
should display a Reynolds number dependence, propor¬ 
tional to 

Both the wavenumber scaling of the temperature spec¬ 
trum and the Reynolds number dependence of the inte¬ 
grated quantities depend thus strongly on the scaling of 
E^(k). The scaling obtained for Eg^\k), 02^^^ and ^4) jg 
identical to the one observed in closure simulations [3. In 
the following it will be investigated by DNS which scaling 
is observed when the discretized Navier-Stokes equations 
are numerically solved. To evaluate the Reynolds num¬ 
ber dependence, we will fix the forcing lengthscale L and 
vary the viscosity. 
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III. SIMULATIONS OF STATISTICALLY 
STATIONARY ISOTROPIC TURBULENCE 
WITH FRICTIONAL HEAT PRODUCTION 


A. Numerical method and required accuracy 


We solve the Navier-Stokes equations with a large-scale 
forcing term, 


dui dui 

dt dxj 


1 dp d^Ui 
pdxi ^ dx'j 


+ fij 


(32) 


and the equation of the total (average plus fluctuating) 
temperature, 


90 90 9^0 V f dui dui dui duj\ 

dt dxi ^ dxf Cp dxj dxj dxi J ’ 


in a cubic three-dimensional periodic domain of size 27r, 
using a standard pseudo-spectral solver. The forcing is 
introduced by a negative viscosity, acting on the modes 
with wavenumbers smaller than 2.5. The initial temper¬ 
ature field is zero and the initial velocity field consists of 
random noise. All results are evaluated once a statisti¬ 
cally steady state is reached were the velocity and tem¬ 
perature fluctuate with a constant variance, and 0 ^ t. 
Once this state is reached, results for (P-^ eg, uf and e 
are obtained by averaging over a sufficiently large time- 
interval to have the errors induced by the finiteness of 
the averaging time to be smaller than errors induced by 
the discretization of the domain. In order to perform this 
time-averaging over a long enough time-interval (70 inte¬ 
gral timescales T = Cfem/e, with Cfcm the kinetic energy), 
the simulations were carried out at a moderate maximum 
resolution of 256^ grid-points. 

It was observed that the results were very sensitive to 
the resolution of the simulations. This issue is well known 
for DNS studies of dissipation-range quantities [mu. 
Whereas the kinetic energy and dissipation were reason¬ 
ably well resolved for simulations with ~ 1, the 

temperature statistics needed about twice this resolution. 
Simulations were therefore carried out using k^a.xV ~ 2. 
The discretization errors were estimated by comparing 
the measured quantities at the highest resolutions at 
fcmax?? ~ 2 and 1.5, respectively, £/ = I /2 - / 1 . 5 I/I/ 2 I 
for a given Reynolds number Rx = 61. The errors in the 
kinetic energy, dissipation rate, scalar variance and scalar 
dissipation were hereby estimated to be 1%, 2%, 14%, 9% 
respectively. These errors will be indicated by error-bars 
in Figure SI Even though these relative errors are quite 
large, the scaling behaviour of the wavenumber-spectra 
was not significantly affected. 


B. Results on the scaling of the velocity and 
dissipation rate spectra 

As was first pointed out by Taylor [l3|, the dissipa¬ 
tion rate becomes independent of the viscosity at high 



FIG. 1: Energy spectra at different Reynolds numbers, nor¬ 
malized using the quantities e, L. 


Reynolds number. The dissipation rate depends then 
only on the large scale quantities L and U where U is the 
root-mean-square value of a velocity component. The 
relation is 


e ~ 


El 


(34) 


where the proportionality constant is independent of the 
Rey nolds number but depends on the large scale flow 
[iq . Since the dissipation rate depends only on the large 
scales, the energy spectrum m, normalized by the quan¬ 
tities e and L, 


= ^ 575^73 ~ ( 35 ) 

should therefore collapse at high Reynolds numbers in 
the forcing- and the inertial-range, if the large-scale forc¬ 
ing mechanism is not a function of the Reynolds number. 
This is shown in Figure [1] to be the case already at the 
relatively low values of the Reynolds number considered 
in the present work. We note that often in literature 
the scaling of spectral energy distributions is displayed in 
Kolmogorov-variables, i.e., length- and timescales based 
on e and v. This allows a collapse of the spectra in the 
dissipation range. However, in the present case, the use of 
the Tayor-variables L and U to collapse the normalized 
spectra in the forcing scales is more efficient to decide 
which scaling describes better the spectra of the dissipa¬ 
tion rate and heat fluctuations. 

The spectrum of the dissipation rate fluctuations is, 
according to expression (ED, proportional to P. If the 
large scales of the velocity field are independent of the 
Reynolds number, as suggested by the results in Figure 
[T](a), this would imply that the normalized spectra would 
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(a) 



(b) 

FIG. 2: (a) Dissipation fluctuation spectra at different 

Reynolds numbers, normalized using the quantities e, L, u, as¬ 
suming expression m to hold, (b) The same spectra nor¬ 
malized according to expression Ill- 

scale as 

E(^\kL) ^ (36) 

and collapse in the large scales. It is shown in Figure 
[H (a) that this is not the case. According to expression 
d, the spectra should rather be normalized as 

^^){kL) = ^^ (37) 

It is shown in Figure [2] (b) that this gives a good col¬ 
lapse of the data. Also, the value of is of order 1/3, 
but the precise value is hard to determine at these val¬ 
ues of the Reynolds number. From these results we can 



(a) 



(b) 

FIG. 3: Temperature spectrum for different Reynolds num¬ 
bers. (a) Normalized using expression (1381) 1. (b) Normalized 
using the expression (EH). 


conclude that the model proposed by Yaglom is suffi¬ 
ciently accurate to roughly describe the Reynolds number 
dependency of the dissipation-rate fluctuation spectrum 
observed in the present work. 


C. Results on the scaling of the temperature 
fluctuations 

In order to verify which of the predictions (I20p or m 
is correct, we use the same arguments as used in the 
previous section for the spectrum E^(k). Indeed, nor- 
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malizing the scaling dini) we should have 

„2 rl/3 

Ei^\kL) = (38) 

Using scaling (EB, we have 

Ej;^\kL) = -^^^Ef(k) ^ {kLY^I^ for a = 0,(39) 
- (fcL)-’5/3+'^ for a = 2/3.(40) 

Clearly, comparing Figures |3] (a) and (b) we see that the 
results are in far better agreement when using the second 
scaling relation. In particular the collapse at small values 
of kL is significantly better using scaling or (gni) 
than for scaling (l38l) . The question whether the inertial 
range timescale is a function of the wavenumber or not 
is not easily answered from this representation. Indeed 
the spectrum does not seem incompatible with a scaling 
proportional to {kL)~^/^'^^ (we have added a line with 
a powerlaw exponent of —1.3), but the inertial range is 
not sufficiently large to give a conclusive answer. 

In FigureUwe show the Reynolds number dependence 
of 0^ and eg, using the normalizations proposed in expres¬ 
sions dsni) and m- The normalized variances should, 

according to our predictions, become independent of the 

~ ( 2 ) 

Reynolds number. This is very clearly the case for , 
which confirms the results observed for the normalized 
spectra. We recall that this quantity is not sensitive to 

the shape of the timescale. For the dissipation of the 

~(2) 

temperature fluctuations, the normalization Cg ' is closer 

to a constant value than e‘'g\ but a clear Reynolds num¬ 
ber dependence is observed. We recall that according to 
the prediction using the integral timescale (a = 2/3) in 
the expression for the heat-fluctuation production-term, 
the normalized temperature-dissipation should scale as 
^3/i/2 ^ ^ powerlaw through the data-points 

in Figure [4] b gives a dependence of correspond¬ 

ing to a value of ^ = 0.37. Note that this value of /r 
gives a scalar spectrum proportional to k~^'^, as indi¬ 
cated in Figure [3] (b). It seems that the non-Gaussian 
model combined with the integral timescale for the La- 
grangian correlation time T{k) is capable of explaining 
all the numerical results. 


IV. DISCUSSION 

The intermittent character of the dissipation rate fluc¬ 
tuations affects dramatically the spectrum of the tem¬ 
perature fluctuations. The wavenumber dependence 
of the temperature spectrum is shown to be strongly 
correlated at large scales, whereas the Gaussian es¬ 
timate and closure expressions of the Eddy-Damped 
Quasi-Normal Markovian type mispredict this spectrum 
to have a correlation-length of the order of the dis¬ 
sipation scale. Our results further indicate that the 


10 '’ 


10 “= 


10 ' 


10 '' 


10 ' 


10 ' 


( 1 ) —^ 
(2) ^-x- 
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10 " 


10 ' 
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10 ' 


10 ' 


(1) —X- 

(2) ^-K- 


10" 


(b) 


FIG. 4: Reynolds number dependence of the normalized tem¬ 
perature variance and the normalized dissipation of temper¬ 
ature fluctuations. Normalizations are given in expressions 
IMl) and (|2^. 


timescale governing the dynamics of the heat-production 
is the integral timescale and not a timescale of the form 
The intermittency coefficient appearing in 
Yaglom’s model for the dissipation-rate fluctuations is of 
the order /i 0.3, if it is estimated from the Reynolds 
number dependence of the temperature dissipation. 

One could speculate on the underlying physical pro¬ 
cesses in terms of flow structures to explain the observed 
results. Indeed, vortices are elongated structures which 
have a diameter roughly between the dissipation length- 
scale and the Taylor lengthscale, depending on the pre¬ 
cise definition, and a length which can extend up to to 
the integral lengthscale. It is tempting to relate this last 
lengthscale to the correlation-length of the temperature 
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(b) 

FIG. 5: Visualizations of (a) the temperature field (red isosur¬ 
faces correspond to positive heat fluctuations: OjOrma = 0.83; 
the blue ones to negative fluctuations QjQrma = —1.26.) 
(b) the vorticity field (iso-surfaces of the enstrophy for 
= 2.77). The Reynolds number is R\ = 77. Visu¬ 
alizations by VAPOR [l^ . 


fluctuations, arguing that they are produced in the strong 
velocity gradient regions around these vortices. 

Another phenomenological explanation could start 
from the observation that clusters of vorticity are cor¬ 
related over the integral lengthscale, a picture which was 
used to model the intermittent character of the pressure 


fluctuations and dissipation-rate fluctuations in [l^ . In 
order to rigorously investigate such mechanisms, precise 
criteria should be introduced to define structures and 
clusters, and this is beyond the scope of the present 
work. However, as a first step, we present in Figure[5]flow 
visualizations, where we show iso-vorticity surfaces and 
surfaces of temperature fluctuations. The main observa¬ 
tion is that, in agreement with the spectral distribution 
of temperature fluctuations, the volume rendering shows 
large-scale correlated smooth temperature iso-surfaces. 
The enstrophy fluctuations show their typical worm-like 
small-scale structure. Comparing figures (a) and (b), a 
correlation seems to exist between regions of intense vor¬ 
ticity clustering and positive temperature fluctuations, 
but more quantitative measures are needed to assess this. 

With respect to the predictions of closure 0, we can 
estimate what the difference in realistic flows will be for 
the temperature variance. Combining expression (IMl) 
with ((Ml) and ((Ml) we find that the variance according to 
the two different model predictions scales as 




(41) 


If we consider a perfectly insulated flow, the mean tem¬ 
perature will, in a steady state, increase linearly in time, 
since 


dt Cp' 


(42) 


The timescale over which the temperature increases, nor¬ 
malized by the temperature variance is then, according 
to ((M() . 


^ t2/3 -1/3 

dO/dt 


(43) 


which corresponds to the integral timescale. This means 
that the rms fluctuations will be of the order of the mean 
temperature-increase over one integral timescale. In the 
closure prediction, the timescale in (1^ was of the order 
of the Kolmogorov-scale, so that for a moderately tur¬ 
bulent flow, the rms fluctuations were expected to be a 
factor Rx smaller in a practical situation and measure¬ 
ment of the temperature fluctuations seemed out of reach 
in a turbulent flow for that case. According to the pre¬ 
dictions of the present work, such measurements will be 
possible, though challenging. 

A question with respect to turbulence theory is how 
these effects could be captured by statistical theory de¬ 
rived from the Navier-Stokes equations. In particu¬ 
lar, which assumptions should be modified in the Di¬ 
rect Interaction Approximation to capture such effects. 
Progress on this seems to have been small since the later 
works of Kraichnan . 
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